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Abstract: 

Earth’s modern atmosphere is highly oxygenated and 1s a remotely detectable signal of 
its surface biosphere. However, the lifespan of oxygen-based biosignatures in Earth’s 
atmosphere remains uncertain, particularly for the distant future. Here we use a 
combined biogeochemistry and climate model to examine the likely timescale of 
oxygen-rich atmospheric conditions on Earth. Using a stochastic approach, we find that 
the mean future lifespan of Earth’s atmosphere, with oxygen levels more than 1% of the 
present atmospheric level, is 1.08 + 0.14 billion years (lo). The model projects that a 
deoxygenation of the atmosphere, with atmospheric O» dropping sharply to levels 
reminiscent of the Archaean Earth, will most probably be triggered before the inception 
of moist greenhouse conditions in Earth’s climate system and before the extensive loss 
of surface water from the atmosphere. We find that future deoxygenation is an 
inevitable consequence of increasing solar fluxes, whereas its precise timing is 
modulated by the exchange flux of reducing power between the mantle and the ocean— 
atmosphere~crust system. Our results suggest that the planetary carbonate-silicate cycle 
will tend to lead to terminally CO,-limited biospheres and rapid atmospheric 
deoxygenation, emphasizing the need for robust atmospheric biosignatures applicable 
to weakly oxygenated and anoxic exoplanet atmospheres and highlighting the potential 
importance of atmospheric organic haze during the terminal stages of planetary 


habitability. 


Introduction 

The field of exoplanet astronomy is on the 
verge of the first detailed characterization of the 
atmospheres of extrasolar planets, with 
heightened expectations for the potential 
detection of spectroscopic indicators of 
habitability, and possibly life, using the next 
generation of space and ground-based 
telescopes!*. However, a robust framework for 
remote life detection requires a detailed 
understanding of the factors that lead to the 
emergence and long-term maintenance of 


atmospheric biosignatures, including a full 


understanding of potential atmospheric 
composition in the broader context of 
planetary biogeochemistry?. Although 


numerous potential atmospheric biosignatures 
have been proposed’, molecular oxygen (Oz) 
and its photochemical by-product ozone (Os) 
remain at the forefront of approaches toward 
remote life detection due to their strong 
absorption at visible and near-infrared 
wavelengths, the importance of biology in 
controlling their atmospheric abundances, and 
the potential for homogeneity in abundance 
throughout the atmospheric column*®. As a 
result, significant effort has been directed 
toward understanding atmospheric 
oxygenation on habitable worlds and the 
factors that might lead to abiotic oxygen-rich 
atmospheres’. 

Earth’s photosynthetic biosphere currently 
supports large gross fluxes of Oz to the ocean- 
atmosphere system—roughly 9X 10!° mol O2 y- 


l1 (ref.!°}—with the result that the modern 
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atmosphere is 20.5% O2 by volume. However, 
the presence of oxygenic phototrophs on its 
own may not be enough to maintain a strongly 
oxygenated atmosphere. Indeed, it is generally 
thought that the abundance of Og in Earth’s 
atmosphere has been well below that of today 
for most of Earth’s history!!, and that the 
modern atmospheric O2 abundance developed 
only after the emergence of land plants!?:!, 
which have evolved to accelerate the 
geochemical cycles of bioessential elements 
(most importantly phosphorus) in Earth surface 
environments!*!5, Thus, charting the history of 
planetary oxygenation after the evolution of 
oxygenic photosynthesis requires an 
understanding of how biospheric productivity 
(and the subsequent burial of reduced materials 
into sediments) changes in response to evolving 
environmental conditions. 

Previous work on the future lifespan of 
Earth’s biosphere has focused principally on 
the links between secular changes in solar 
luminosity, the stability of the carbonate- 
silicate geochemical cycle, and the loss of 
surface water to space!®?9, One robust result 
that emerges from this theoretical framework is 
a continuous decrease in the abundance of 
atmospheric CO2 through time, as steadily 
increasing solar luminosity inexorably drives 
down the abundance of atmospheric CO» that 
is required to balance rates of crustal 
weathering with volcanic CO: output!®!77!, 
However, a number of models suggest that in 
spite of this thermal buffering mechanism 
Earth will enter the “moist greenhouse” climate 


state at some point within the next ~2 billion 


years (Ga)—a regime in which elevated 


atmospheric temperatures allow  Earth’s 
stratosphere to become moist, leading to 
permanent water loss from Earth’s surface!™ 
19,22,23. In this 


atmospheric CO2 may result in a fundamental 


addition, drawdown of 
shift to a photosynthetic biosphere that is CO2- 
limited—for instance, some models predict that 
terrestrial C3 plants will cease to be viable at 
Earth’s surface less than ~500 Myr into Earth’s 
future!®. If true, this might place a long-term 
physical limit on the ability of photosynthetic 
maintain high levels of 
giving 

trade-off between 


biospheres to 


atmospheric oxygen, rise to a 


fundamental long-term 
stellar evolution, the geologic carbon cycle, and 
the intrinsic timescale of atmospheric 
oxygenation. 

Here, we examine the future lifespan of 
Earth’s oxygenated atmosphere using an Earth 
system model of biogeochemistry and climate 
(Fig. 1; see Methods). Assessing the lifespan of 
Earth’s 


quantitative models that include a complete 


oxygenated atmosphere requires 
treatment of oxygen-related biogeochemical 
processes at the planetary surface over geologic 
time, as well as mechanistic linkages between 
surface environments and the planetary 
interior. Our model builds upon previous 
similar Earth system models!? that track the 
coupled carbon, oxygen, phosphorus, and 
sulphur cycles in the exogenic (i.e., atmosphere, 
ocean, and crust) system (Fig. 1). However, we 
extend this framework by incorporating a 
global biogeochemical methane (CH4) cycle, 
of key biological 


including a range 
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parameterized atmospheric 


Oo-O3-CH4 


system?#?9, and the radiative impact of CH4 on 


metabolisms, 


photochemistry within the 


global energy balance?®. 
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Fig. 1. Schematic model structure. Boxes denote 
reservoirs, whereas arrows denote flux terms. The 
model tracks the major reservoirs and transfer 
fluxes within the surface carbon (C), sulphur (S), 
oxygen (O), and phosphorus (P) cycles, along with 
a comprehensive treatment of ocean 
biogeochemistry, and long-term transfers between 
the crust-ocean-atmosphere system and the mantle. 
DOA = degree of anoxia. 


Our model is also the first of its type to 
explicitly evaluate redox exchange between the 
crust and mantle, allowing a comprehensive 
evaluation of the planetary processes 
controlling atmospheric Oz levels over geologic 
time. The model is designed to capture the 
Earth 


biogeochemistry and climate, but is abstracted 


major components of system 
enough to allow for large model ensembles that 
can be run for billions of years of model time. 
We employ a stochastic approach in which we 
draw randomly from a wide range of potential 
values for parameters controlling the major 


geophysical boundary conditions of the model, 


including variability (cycle and amplitude) in 
outgassing rates and erosional forcing (see 
Methods and Extended Data Figs. 1-3), and 
initialize the model at 600 Ma. We then 
repeatedly run the model forward in time until 
the present day (n =~400,000), and subsample 
model runs (n = 4,787) that roughly reproduce 
modern atmospheric composition (pO. = 
21+1% and pCO:2 = 3004150 ppmv), global 
average surface temperature (288+2K), and 
seawater sulphate (SO4*) levels (SO47 = 25415 
mM). The simulations from this subsampled 
ensemble are then run forward in time in order 
to estimate the future lifespan of Earth’s 


oxygenated atmosphere. 


Solar brightening promotes atmospheric 


deoxygenation 


Our default analysis (Fig. 2 and Extended 
Data Figs. 4-8) indicates that atmospheric O2 
levels will decrease substantially in Earth’s long- 
term future (Fig. 2a). Whilst the magnitude 
and timing of this trajectory show a certain 
degree of uncertainty, the model-projected 
decrease in atmospheric Og abundance is 
robust and suggests that the future lifespan of 
Earth’s oxygenated atmosphere will be less 
than 1.5 Gyr. Specifically, when we define the 
future lifespan as the period within which 
atmospheric O2 remains above a given 
threshold value, the lifespan of atmospheric O2 
than 1% 
atmospheric level (PAL)—11%,paL—1s estimated 
at 1.084£0.14 Gyr (10). Increasing the threshold 
value to 10% PAL results in a broadly similar 


lifespan of T 1o~paL = 1.0540.16 Gyr (lo), a 


levels greater of the present 
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consequence of the extremely short timescale 
over which atmospheric deoxygenation occurs. 
Several factors are responsible for the secular 
decrease in atmospheric O2 (see 
Supplementary Information), most notably 
increasing surface temperatures and 
atmospheric CO» starvation of the terrestrial 
photosynthetic biosphere driven by solar 
brightening (Fig. 2c and Extended Data Figs. 
4-7). Indeed, when we repeat our analysis with 
a constant solar luminosity we find that no 
secular O2 trend remains (Fig. 2a and 
Extended Data Fig. 4). As the geologic drivers 
of the decreasing COs trend (solar brightening 
and mantle cooling) are largely undisputed, the 
first-order features of our projected 
atmospheric O» trajectory should be robust. 
Our model further predicts that CO2 
limitation of Earth’s photosynthetic biosphere 
should lead to photochemical 
destabilization of Earth’s 


atmosphere and an abrupt shift in atmospheric 


ultimately 


oxygenated 


Oy abundance to very low values (Fig. 2a). 
Indeed, our model projects an atmospheric 
composition for Earth’s long-term future that is 
similar in many respects to that of the Archaean 
Earth prior to the so-called “Great Oxidation 
Event’?’. In particular, atmospheric O2 drops 
to values that are many orders of magnitude 
below that of the modern atmosphere (<10-6 
PAL), while the atmospheric abundance of 
CH; increases sharply (Fig. 2b). However, a 


critical difference from the Archaean Earth 
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Fig. 2. Evolution of atmospheric chemistry in our stochastic analysis. Blue, orange, and green lines 
represent atmospheric O2 (a), CH4 (b), and CO2 (c) concentrations obtained from our default model ensemble 
(n = 4,787), along with the predicted atmospheric CH4/COz ratios (red lines in d). Horizontal dotted lines in d 
represent a range of values between 0.1-0.2, broadly considered to define a threshold for the inception of organic 
haze formation?**°. All panels also show results obtained by running the same model ensemble with constant 
solar luminosity (grey lines). Phan = Phanerozoic. 


system is the prediction of much lower 
atmospheric CO% levels in Earth’s long-term 
future (Fig. 2c) and thus the potential for 
greatly elevated CHs/CQOz ratios (Fig. 2d). 
Such elevated CH4/COz ratios are generally 
predicted to result in the formation of ‘organic 
haze’ in the atmosphere?®“°, with potentially 
important impacts on climate, atmospheric 
redox chemistry, and remote life detection. 
Our results thus suggest that organic haze will 
be a critical component of Earth’s long-term 
climate stability into the far future, as well as a 
potentially promising atmospheric biosignature 
on Earth-like worlds round late main-sequence 
stars such as the exoplanet candidate Kepler- 
452b (orbiting a G2 star with an age of ~6+2 
Gyr and receiving ~10°% more insolation than 


Earth receives from the Sun today)*!°?. 


Controls on the timing of atmospheric 


deoxygenation 


Our default analysis provides the first 
systematic estimate of the potential lifetime of 


Earth’s oxygenated atmosphere into the far 


future. However, it is important to note that our 
default model ensemble includes a robust 
terrestrial biosphere, which may not be a 
universal feature of habitable Earth-like planets. 
To examine the impact of the terrestrial 
biosphere on the longevity of the O2/O3-based 
biosignatures, we performed an additional 
ensemble analysis in which we remove the 
terrestrial biosphere from our default successful 
runs (Extended Data Figs. 4 and 10). As 
expected, the absence of land plants results in 
somewhat lower atmospheric Oz levels 
throughout the course of planetary evolution 
relative to our default analysis. However, we 
find that the future lifespan of Earth’s remotely- 
detectable oxygenated atmosphere is similar to 
that of our default analysis (e.g., Ti%PAL = 
1.04*0-17 0.16 Gyr; cf. 1.08+0.14 Gyr for the 
default simulations) (Fig. 3a). This is an 
important result, as it strongly suggests that the 
presence or absence of terrestrial biosphere 
exerts a secondary control on the timescale of 
the planetary deoxygenation and long-term 


atmospheric redox state. Although it is difficult 


to explicitly model uncertainties associated with 


biological evolution and adaptation in 


photosynthetic lineages, our statistical 
approach explores a wide range of parameter 
space for key biological response functions (see 
Methods and Extended Data Figs. 1-3), with 
results suggesting that uncertainty in the 
parameters governing the response of the 
biosphere plays a secondary role in controlling 
the future lifespan of oxygenated atmospheric 
conditions on Earth (Extended Data Fig. 9). 
Our analysis also delineates an important 
relationship between the future lifespan of 
oxygenated atmospheric conditions and the 
long-term average flux of reducing power from 
o,(Red) 
(Fig. 3b). In particular, our model suggests that 


the mantle to the surface system, &® 


larger @, (Red) values tend to yield shorter 


oxygenated lifespans. Improvements in 


quantitative estimates of subduction and 
mantle degassing over time would reduce 
uncertainty in estimates of the future lifespan of 
Earth’s oxygenated atmosphere, and will be 
critical for extending our model framework to 
very different tectonic regimes. In addition, our 
analysis suggests that oceanic redox 
chemistry—which ultimately controls oceanic 
P levels and the redox state and subduction 
fluxes of major C, S, and Fe species in the 
exogenic system—would be critical for the 
long-term trajectory of atmospheric redox state 
on habitable Earth-like planets (Fig. 3b and 
Extended Data Fig. 4). The lifespans of oxygen- 
rich atmospheric states on Earth-like planets 


are thus likely to be modulated by mantle redox 
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Fig. 3. The future lifespan of Earth’s 
oxygenated atmosphere. (a) Posterior 
probability density distributions from our default 
analysis for Earth systems with (dark grey) and 
without (light grey) a terrestrial biosphere. Open 
circles denote median values, while the error 
bar and shaded region denote 1o and 95% 
credible intervals, respectively. (b) The future 
lifespan of Earth's oxygenated atmosphere 
(t1%PaL) as a function of the average net input 
flux of reducing power from the mantle to the 
surface system. Positive values indicate net 
input of reducing power from the mantle to the 
surface system (see Methods). Each point is a 
single model run, while the colour scale denotes 
data density within the overall model ensemble. 
We note that there are a small number of 
scenarios in which atmospheric O2 levels do not 
reach 1%PAL during the simulation, but such 
scenarios are very rare (~4%) and would rely on 
a fortuitous combination of parameter values. 


state, rates of surface-interior exchange, and 
their evolution through time. 

It is important to note that there are 
multiple biogeochemical and climate processes 
that are not considered in our model that may 


play a role in constraining the future lifespan of 


Earth’s biosphere and the timing/mode of 
transition to more reducing atmospheric 
conditions. In particular, “reverse weathering” 
(the formation of authigenic silicates in marine 
sediments, resulting in net COz release to the 
ocean-atmosphere system)?* could potentially 
extend the lifespan of oxygenated atmospheric 
conditions under certain scenarios by 
which 


atmospheric COs is above the levels expected 


prolonging the timescale over 
to result in COs limitation of the photosynthetic 
biosphere*#. On the other hand, this process 
global 


temperature in the face of increasing solar 


would commensurately elevate 
luminosity, with a corresponding potential for 
limitation of biospheric Oy fluxes (Extended 
Data Fig. 1). The net effect on the future 
lifespan of Earth’s oxygenated atmosphere is 
uncertain, but is likely of secondary importance 
to the overall trend presented here (see also 
Extended Data Fig. 5 for the sensitivity to 
terrestrial weatherability). In addition, we 
hypothesize that haze-induced climate 
cooling?! could potentially act as a brake on the 
overall magnitude of atmospheric 
deoxygenation, or result in the inception of 
oxygenation/deoxygenation cycles during 
Earth’s terminal habitability. Both of these 
topics represent important areas of future 


research. 


Implications for the search for life on 


exoplanets 


Our 


implications for the search for life on Earth-like 


results also have important 


planets beyond our solar system (e.g., habitable 
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planets with abundant liquid water at the 
surface, exposed silicate crust, and a biosphere 
with oxygenic photosynthesis). From the 
perspective of planetary evolution, our results 
imply that atmospheric oxygenation is not a 
permanent condition on habitable worlds 
hosting oxygenic photosynthesis, and that only 
of Earth’s 


characterized by robustly detectable levels of 


a fraction history will be 
atmospheric O2 (Fig. 4). For example, it is 
estimated that a G-dwarf stellar flux ~20% 
greater than that of the present Sun (S/So ~ 
1.2) is required to achieve a water loss timescale 
of ~1 Gyr in the moist greenhouse regime”. 
Assuming it will take ~2 Gyr to achieve S/So ~ 
1.2 (ref.°), together with the assumption of 
surface liquid water by 4.4 Ga, yields an 
estimated habitable lifetime on Earth of ~7.4 
Gyr. Our results indicate that direct detection 
of Os 
challenging for all but ~1.5-2.0 Gyr of this 
history — or roughly 20-30% of Earth’s 


at visible wavelengths would be 


lifetime as a habitable world*®. On the other 
hand, observations of Os at UV wavelengths?» 
could potentially significantly extend this 
observability timescale (Fig. 4). In any case, our 
model suggests that a fundamental feature of 
the carbonate-silicate cycle, considered central 
to the long-term maintenance of habitable 
climate conditions on terrestrial planets, may 
ultimately act to drive atmospheric 
composition toward an anoxic, hazy state 
similar to that of the Archaean Earth but with 
CO2. These 


predictions emphasize the need to better 


much lower atmospheric 


understand the factors regulating redox 


limitation 
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Fig. 4. The coupled evolution of Earth’s 
biosphere and atmospheric chemistry. 
(a) Solar constant (normalized by modern 
value)'®. (b) Limiting factor for biospheric 
activity. (c) Global net primary production 
(normalized to the modern value)***°. (d) 
Atmospheric CO2 level (in bar) *’. (e) 
Atmospheric O2 level (relative to the present 
atmospheric level, PAL)‘'. (f) Atmospheric 
CHa level (in parts per million by volume, 
ppmv). The atmospheric CH4 level and 
global NPP during the mid-Archean and 
mid-Proterozoic are based on ref.4°4°. The 
historical trajectory of biospheric productivity 
suggests that despite evolving limitations on 
overall biospheric fertility the biosphere has 


C4 plants 


significantly altered the chemical 


-4 3 -2 -1 0 1 
Time (Gyr) 


balance between the planetary interior and 


surface environmentst®-44, and provide 
additional motivation for implementing a wide 
wavelength range in future exoplanet 
characterization missions as a contingency 
against “false negatives” for remote life 
detection?:38, 

In sum, our stochastic analysis suggests that 
the of Earth’s 


atmosphere is a robust outcome of increasing 


eventual deoxygenation 


solar luminosity, irrespective of large 


uncertainties in geophysical/biological 
boundary conditions, and that the collapse of 
Earth’s oxygenated atmosphere is likely to 
occur before the inception of moist greenhouse 


conditions in Earth’s climate system although 


composition of the atmosphere*®46485°, Our 
model, however, demonstrates that Earth’s 
modern biosphere reaches the turning point 
at which the activity level of the biosphere 
begins to decline (see also Extended Data 
Figs. 4-7), and that the availability of 
inorganic carbon in combination with 
increased surface temperature will become 
the major limiting factor for global 
ecosystems. The orange-red shaded region 
> denotes the inception of the moist 
greenhouse climate regime. 


its precise timing is modulated by net input flux 
of reducing power from the mantle to Earth’s 
In addition, 


projections of atmospheric chemistry implicate 


surface system. our model 
organic haze as a critical factor in regulating 
Earth’s long-term future climate. Evaluating 
these factors using more sophisticated 
biogeochemistry and climate models is an 
important task for future research. 
Nevertheless, our analysis clearly identifies the 
first-order response of atmospheric O2 to 
declining CO2 and biospheric productivity, 
opening a new perspective on the future 
evolution of Earth’s climate and atmospheric 


redox chemistry. 
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Methods 


Model description. Our model is based on the Carbon Oxygen Phosphorus Sulphur Evolution 
(COPSE) model!?:*8, the default version of which is described in full elsewhere!?:*8. The model has 
been extensively tested and validated against geologic records from the most recent ~500 million 
years of Earth’s history and can predict the coupled histories and controls on atmospheric COs, 
Oz and ocean composition (PO4-, NOs, SO4?) on geologic timescales. The COPSE model was 
recently revised by ref.°!, but the version used here is based mainly on a previous version of the 
model described in ref.!?!5, because we have constructed the model with the aim of keeping the 
design as simple as possible. We do not anticipate our overarching conclusions to be altered by the 
differences. Our aim here is to produce a quantitative model that comprehensively encapsulates a 
range of key biogeochemical processes in oxic/anoxic biospheres and makes resulting predictions 
for the future lifespan of Earth’s oxygenated atmosphere. Specifically, we made the following 
additions to the basic structure of the model: (1) metabolic CH4 production/ consumption processes 
in the ocean, (2) parameterized atmospheric photochemistry within the O2-O3-CH4 system?**», (3) 
the radiative impact of CH4 on global energy balance?®, (4) a dependence of marine phytoplankton 
growth rate on surface temperature?! and marine inorganic carbon availability, (5) exchange of C, 
S, and Fe between the exogenic (ocean-atmosphere-crust) system and the mantle, and (6) a full 
redox (O2) budget in the exogenic system. 


Primary production and decomposition. In the COPSE model, the activity levels of terrestrial and 
marine ecosystems are evaluated separately. The global net primary production (NPP), Jnpp (in 
terms of organic carbon), is given by: 


— yocn land 
J pp ~ J pp + J pp, (1) 


where Jnpp°® and Jnppl#4 denote marine and terrestrial biospheric productivity, respectively. 
Jnpp°™ is assessed as follows: 


Ixpp = Sy Se J Npe., (2) 
where * denotes the modern value of ~45 Gt C yr7!, or 3750 Tmol C yr! (ref.!%), and fh and from 


denote the effect of nutrient availability and surface temperature on oceanic productivity, 
respectively. A Liebig’s minimum law is adopted for fn: 


Minimum ren , Top P, M (3) 


* 
Ton N 


CO, +HCO; | 


> 


i= 


where N, P, and Mco2+Hco3- denote reservoir sizes of seawater nitrate, phosphate and inorganic 
carbon (CO2aq and HCOs’) available for phytoplankton growth, and rc,p (= 117) and rns (= 16) 
denote the “Redfield ratios” of photosynthetic biomass. In the original COPSE model the growth 
rate of marine phytoplankton is limited by the availability of nutrients (N or P). However, even in 
the modern ocean, CO+-limitation can reduce growth rates among some phytoplankton taxa°??. 
Given a secular decrease in atmospheric CO» levels (Fig. 2c), marine phytoplankton could be 
limited by the availability of dissolved inorganic carbon (CO2aq and HCQOs3;) in the remote future. 
The dissolved concentrations of CO2aq and HCO3° are calculated in the model by solving the 
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seawater carbonate system (see Supplementary Table 3). For the temperature response of marine 
phytoplankton growth, we adopt the following general function based on ref.*4; 


0 (T, < 273 K) 
on (TP) =4 kd] UEA (273 K <T, <301 K) 
T s27 TA Seg (4) 
a jy OE 
k| en (cae | a L (7, 2301 K) 
AT pt Tot 


where k denotes a normalizing factor (= 1.2748), and Top: and ATopt are assumed to be 301 K and 
301 K3, respectively. The maximum temperature for phytoplankton growth is controlled by Tres, 
whose uncertainty is explored by the statistical analysis (see below and Extended Data Figs. 1 and 
3). 


Almost all organic matter produced in the surface ocean is decomposed before burial. Fractions of 
organic matter decomposed by aerobic respiration (combined with denitrification/nitrification and 
microbial sulphate reduction/sulphide oxidation) and methanogenesis are expressed as follows: 


O 
fo, = O+Ko,: (5) 
Jon, =1- fo, (6) 


where O denotes the oxygen reservoir size in the combined ocean-atmosphere system and Ko2 = 
13.6X10!8 mol. The total organic carbon decomposition rate is given by pow, — phon = pron, 


org org 
where Jorg’°™ denotes the burial rate of organic carbon in marine sediments that is proportional to 


the square of the normalized primary production as in the original COPSE model**®. The global 
methane generation rate via methanogenesis ( 2CH,O —>CH + CO,) is given by 0.5x fon * IRS. A 
fraction of this methane is consumed by aerobic/anaerobic methane oxidation processes (aerobic 
methanotrophy and AOM), given by: 


O 


7 O+Ko, ) 


where Kog = 0.273X10!8 mol (ref. 54). Given the stoichiometric relationships, the net production 
rate of CH4, COs, and O» in the ocean interior are given by: 
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hn _ =3(I- ô) Tns Ire (8) 
en s-( fo,+ (148) fr, Jer 0) 


Jo, = NPP Us +5 fox, J Jee”. (10) 


For simplicity, the uncertainty associated with partitioning of CH4 consumption between aerobic 
methanotrophy and AOM is accounted for implicitly by prescribing Koo’. 


The NPP of Earth’s terrestrial biosphere is scaled by a factor V which represents global land 
biomass normalized by the modern value**: 


Tis =v age a1) 


where the NPP of the present terrestrial biosphere, Jnpp®™t*, is ~60 Gt Cyr!, or 5000 Tmol C yr! 
(ref.!°), and V is calculated, as follows: 


V = foy (03): fe ` Frire(O2) Vapp(O2,CO2,7). (12) 


where Vnpp denotes the effects of atmospheric O2 and CO» concentrations and surface temperature 
on terrestrial biomass as in the original model (Extended Data Fig. la,b)*®, fire denotes the effect 
of the fire feedback!?:*®, and fp denotes a forcing factor representing selective-P weathering by 
vascular plant activity!?:!®. Our model introduces an additional factor, fuv, in order to capture the 
effect of UV on the terrestrial biosphere as a function of atmospheric Oz levels: 


(13) 


fa O:)= am| POLPAD 


Cyy 


where pO2(PAL) denotes atmospheric Oz levels in terms of PAL, and cuv is treated as one of the 
model parameters explored during our stochastic analysis (see below). The minimum CO: 
concentration for the terrestrial plant activity is also treated as one of the uncertainties, because the 
ecological COs: threshold is not precisely known and its absolute value is likely to change with 
global mean temperature via its influence on photosynthesizer respiration rate°»®. 


The CH; flux from the terrestrial biosphere is assumed to be proportional to the burial rate of 
organic matter on land, as follows: 


te vie J> land / J® A land 
= org 


land,* 
Se (14) 


where frd denotes the temperature effect on the CH; flux: 
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fr = acu, exp (bon, T/T) (15) 


where To denotes the surface temperature in terms of °C. On the other hand, the net production 
rate of CH4, COs, and O?» in the terrestrial biosphere can be written, as follows: 


an 1 r,lan 
Wi aT 06 
1 
= So, +5 (I+6) son, J, (17) 
~ Te -(go, tgan Je (18) 


where Joland (= Jypp™d — JorgP!and) denotes the decomposition rate of organic matter in the terrestrial 
biosphere, and the fraction of methane oxidation is calculated by equation (7) above. A fraction of 
organic matter decomposed by methanogenesis, gcn4, is determined from equations (14) and (16). 
Then goz is determined from g o, =l- Zon, 


Photochemistry. The revised model includes parameterized photochemistry that allows 
calculation of the concentrations of atmospheric O2 and CH4. We incorporate a scheme for 
parameterized O2-03-CH4 photochemistry based on the previous study of a 1-D photochemical 
model?4, Specifically, the CH4 oxidation rate kcH4ox (in mob! y!) is expressed as a polynomial 
function of the reservoir sizes of O2 and CH; (O and M)”: 


—p/ j J 2 f 3 J 4 J 5 J 6 
log Ker ox =Q) +A Po, tA Po, +3 Po, +4 °Po, tA °Po, TAs Po, , (19) 


where aj is fitting coefficients for a given atmospheric CH; levels and @oz is given by logpOx(bar) 
(Supplementary Table 4). The oxidation rate is evaluated based on Fig. 3 of ref. 24 showing the 
oxidation rate as a function of pOz and pCH4. We took the relationship between kcH4ox and pOe 
for pCH4 of 10-6, 10°, 104, 103, 2x10- bar, and kcH4ox is calculated as a function of pO» and 
pCH4 with a log-linear interpolation method. 


For the atmospheres examined in this study, CH4 and H2O are the dominant hydrogen-bearing 
species in the stratosphere, allowing the hydrogen escape flux, Jesc, to be calculated as: 


J ese z Se à Jr (H,) = Jia F Jess, (20) 


where the proportional coefficient, fesc, is 6.68 10!° mol yr"! and fr(H2) denotes the total hydrogen 
molecule mixing ratio: 


f,(H,) = fH,0+2/CH, (21) 
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Therefore, fluxes of hydrogen escape via CH4 and H2O are written in terms of the abundances in 
the stratosphere as: 


Jie = des ` JH,O 


(22) 
i = 2 : Tee ` JCH, 


We fit the empirical relationship between the stratospheric water vapour concentration and global 
average surface temperature obtained from a suite of 3-D GCM simulations of ref. 3 with a 3*4 
order polynomial in surface temperature: 


log, fH,O= f + fT + fil + fT’: (23) 


where fi denotes the fitting coefficients (fo = 198.34, fi = -1.6856, fp = 4.147x10°%, and f3 = - 
2.7947x10-°) (Supplementary Fig. 2). 


Subduction flux. Despite its crucial role in controlling the redox budget in the exogenic system, 
a process not typically represented in large-scale biogeochemical models is the subduction of redox 
sensitive elements (C, S, and Fe species) into the mantle. The exchange flux between Earth surface 
and the mantle influences not only the sizes of the surficial reservoirs but the oxygen balance of the 
atmosphere on the long timescales considered here. For the modern Earth the subduction flux of 
organic C is estimated at <0.2-0.7 Tmol Oz equiv. yr! (ref. 40-42), and that of sedimentary pyrite 
would also be small (<1 Tmol S yr!)#?°8, On the other hand, the rate of subduction of ferric iron 
could be 3 Tmol Fe yr! (0.75 Tmol O2 equiv. yr-!)*°, or as high as 7 Tmol Fe yr! (1.75 Tmol O2 
equiv. yr!)°9 to 19 Tmol Fe yr! (4.8 Tmol O? equiv. yr!)**. However, there is a large uncertainty 
in the major oxidants (O2 vs. SO4?) for Fe?* during aqueous oxidation of seafloor®?. These redox- 
dependent subduction fluxes probably vary in intensity in response to changes in tectonic activity 
and oceanic redox state. In this study, the subduction fluxes of C, S, and Fe are parameterized as 
follows: 


* J. 
Jog = Jag : a | " Io, (24) 
org 
s s,* Ja 
J = S oxo ` ss } Ja (25) 
carb 
e h 
Joy = Soy [Bs (26) 
py 
* 1- DOA 
Taal P 
Fe-ox Fe-ox (a ) Je ( ) 


where * denotes a reference value, which is treated as a free parameter. The rate of subduction is 
scaled to the tectonic activity represented by a forcing fc (D in the original COPSE model). The 
subduction fluxes of organic carbon and carbonate carbon are also proportional to their burial 
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fluxes. We also assume that the subduction flux of Fe-oxides is a function of the degree of anoxia 
(DOA) which represents oceanic anoxicity. The rate of subduction of ferric iron is thought to be 
sensitive to the redox landscape of the global ocean (e.g., sulfidic or Fe-rich anoxia), but this is 
outside the scope of the present work. Even if subduction of iron oxides is assumed to continue in 
an anoxic, Fe-rich ocean, this would not change our overarching conclusion that the future lifespan 
of the oxygenated atmosphere is less than 1.5 Gyr because the subduction of ferric iron represents 
the loss of oxidizing power from the surface system. Given the large uncertainty in the magnitude 
of subduction fluxes, we explore a wide range of reference values in our stochastic analysis 


(Extended Data Fig. 3). 


Global O2 budget. In this study, redox balance in the model is tracked in units scaled to the 
oxidizing power of O2. The model explicitly accounts for fluxes of oxidants and reductants into 
and out of the surface system (Supplementary Fig. 3). The global redox budget (GRB) for the 
ocean-atmosphere system (OAS) is written as follows: 


GRBoxg = 0.54, 


nose + (Jee +I Te, I) +(2J, -2" )— 2h. J 0.255 


org r Fe-ox 


(28) 


where JHesc denotes the irreversible hydrogen escape to space, and second and third terms represent 
the influence of organic C and pyrite S subcycles on the GRB. Jaos! and Jreat are the mantle 
degassing of H2S and other reducing gases (e.g., Hə). Jre-ox’ is the subduction flux of ferrous iron. 
In this study, we also revised the baseline COPSE model to capture the O2 budget in the exogenic 
(ocean-atmosphere-crust) system on geological timescales, according to: 


GRB, =0.5J 


Hesc 


tJ, te) 2d, 


ea (29) 
The primary sink of reducing power from the exogenic system is the irreversible escape of hydrogen 
to space (JHesc) and the subduction of organic C and pyrite S to the mantle (Jorg and Jpy’), while the 
primary source of reducing power is outgassing of reduced gases from the Earth’s interior (Ju2s* 
and Jre). The subduction of ferric iron (Jre-ox’) is also a loss of oxidizing power (i.e., source of 
reducing power). The net input flux of reducing power from the mantle to the exogenic system is 
defined as: 


O(Red) =—J,,, 25, + 0.25 Tino + iis + Trea, (30) 


C-OX 


In Fig. 3b of the main text, we assess the impact of long-term ®(Red) on the future lifespan of 
oxygenated atmosphere. To do this, we average ® (Red) over the future lifespan: 


I i @(Red)dt 


P, (Red) = 61) 


where ti denotes the future lifespan for the different values of O2 threshold i (in PAL). Existing 

empirical estimates of the mantle degassing of reducing gases via submarine volcanism, 0.70.3 

Tmol O2 equiv. yr’! (ref. 60), the limited O2 consumption through seafloor Fe** oxidation 0.14 

Tmol Oz equiv. yr! (ref. 60), and the subduction flux of organic C of <0.2-0.7 Tmol O2 equiv. yr 
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1, yield a present net input flux of reducing power from the mantle to the surface system, (Red)", 
of >-0.16 Tmol O2 equiv. yr!. A resampling analysis of ®(Red)* from our successful model 
simulations predicts a slightly lower value of -0.24+0.27 Tmol O2 equiv. yr! (Supplementary Fig. 
4). 


Climate model. Global mean surface temperature T, (in K) is estimated with a zero-dimensional 
energy balance model (0-D EBM) in which the absorption of solar radiation is balanced by the 
heat loss by outgoing longwave (1.e., infrared) radiation (OLR) to space (Jorr): 


S 
(1-a):=2 = Jorr RF, (32) 


where @ and § denote the top-of-atmosphere (TOA) albedo and the value of the solar constant 
at the time under consideration, respectively. For Jorr and a, we adopted the parametric 
expressions proposed by ref.°!, with the following variables: p = pCO» (in bar); ọ = In(p/3.3x10- 
+); p = cos(z), where z is the zenith angle; and surface albedo as. a; is written as a function of Ts: 


A (DST) 
C 
a, = Cice (ai Oyo )* (T. -T y (Tice < T; Š D) (33) 
Bic (T, < Tice) 


where it is assumed that Earth is globally ice-covered for Ts < Tice (= 263 K) with an albedo ice 
(= 0.65), and that Earth is characterized by a relatively low albedo ao (= 0.35) for T, > To (= 273K). 
For Tice < T, < To, the albedo is interpolated with a quadratic interpolation method®?. RFcus 
represents the radiative forcing by methane, which is quantified based on ref.. Recent 1-D 
modelling results suggest that relatively small increases in solar luminosity could conceivably trigger 
a runaway greenhouse on Earth®:6*, However, 3-D climate models indicate that Earth’s climate is 
likely to be somewhat more stable than this as a consequence of expansion of unsaturated regions 
in the upper atmosphere attendant to changes in Hadley circulation that occur at higher solar 
fluxes?3-6, Given the quantitative consistency between the above 0-D EBM and recent 3-D climate 
models” (Supplementary Fig. 1), we consider our utilization of 0-D EBM a defensible 
approximation of planetary climate for implementation on 10°-10°-year timescales. Nevertheless, 
it is important to point out that future projections of Earth system behavior will likely be modulated 
by uncertain cloud feedbacks and continental configurations. As a result, in our MC analysis we 
explore uncertainties in OLR and TOA (Extended Data Fig. 3). 


Stochastic analysis. The computational efficiency of the model enables us to examine the future 
lifespan of the oxygenated atmosphere via a stochastic analysis, in which we randomly sample from 
a prior distribution for poorly constrained parameters and quantify whole-model uncertainty 
(Extended Data Fig. 3). We draw poorly constrained initial values for crustal C and S reservoir 
sizes at 600 Ma (million years ago) over a wide range of values. We also vary the reference values 


19 


In press at Nature Geoscience 


of fluxes between the exogenic system and the mantle (subduction and mantle degassing) to capture 
uncertainties in model parameterisation. After running the model forwards in time until the present, 
we only sample model runs that allow for acceptable atmospheric composition (pO2 = 2141% and 
pCO» = 3004150 ppmv), global surface temperature (28842K) and seawater sulphate (S04 
levels (SO4- = 25415 mM). The model runs that fail to obtain reasonable values have been 
abandoned. For our default stochastic analysis, we performed nearly 400,000 simulations to span 
the whole parameter space, and the resultant successful runs were n = 4,787 which were continued 
to be run for estimating the future lifespan of the oxygenated atmosphere. 


For the future evolution of the degassing factor (fc) and erosion factor (fr; U in the original model), 
we assumed following functions: 


je bean 8 , C 64) 


To 


t-T t \" 
=| 1+a, sinl 2 R LII 1+—— (35) 
aR F (2z r | a) 


where the first term on the right-hand side denotes the periodic changes with an amplitude a and 
a period of t (in Myr), and the second term represents the secular decrease in the activity of 
tectonics®®, These relationships have been adopted for examining to the geological past, and we 
consider that this assumption is a reasonable first-order approximation as there is no a priori reason 
for thinking that this style of decay will not continue into Earth’s far future. Nevertheless, we also 
emphasize that this is very under-studied and it is possible that non-linearity in tectonic regime will 
be important to consider for Earth and other habitable rocky planets in future work. 


We draw model parameters from wide distributions that represent the uncertainty in their values 
(Extended Data Fig. 3). For example, the uncertainty in the biological responses is assessed by 
changing biological parameters, such as FERT. The strength of terrestrial weathering feedback is 
also a topic of debate®”-69, and we explore its uncertainty by changing the apparent activation 
energy (ACT) and a runoff factor (RUNsil)’°. The amplification of terrestrial weatherability by 
land plants is also an important uncertainty’!:’*, which we explore here by changing LIFE”. To 
account for the uncertainty in the value of cuv, we draw it from a log-uniform distribution between 
10-25 PAL and 10-15 PAL (Extended Data Fig. 1c). At the higher end of this range, terrestrial plant 
activity is suppressed when atmospheric O2 is <10% PAL. At the lower end, the threshold is 
lowered to ~1%PAL. We consider these two end-member cases to cover the full range of 
uncertainty of UV impacts on the terrestrial ecosystems. We also explore uncertainty in the 
strength of redox-dependent P recycling in marine sediments by changing Corg/Porg ratios for 
anoxic sediments. 


The sampling procedure is run for long enough to obtain stationary distributions for posterior 
probability densities. However, it is possible that our model predictions may underestimate the 
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inherent variability in the future Earth’s evolution. For instance, the possibility of random 
catastrophe, such as the asteroid impacts, could also introduce additional uncertainty in the 
estimate of future lifespan of Earth’s oxygenated atmosphere. We also note that the uncertainty in 
the longevity of human impacts on Earth’s future environment is beyond the scope of this study. 


Data availability. The numerical data for Figure 3 in the main text and Extended Data Figures 
4 and 5 are presented in the Source data. The big data obtained by the statistical analysis is 
available to download at https://doi.org/10.6084/m9.figshare.13487487.v1. 


Code availability. Our Fortran source code is available at https://github.com/kazumi- 
ozaki/lifespan. 
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